###############################################################################################
# 
# Sensitivity Analysis: No Non-Disclosure
# 
###############################################################################################

#########################################################################################
### Democratic Republic of Congo

m.liar <- ictreg(Y ~ 1, treat="treatment", J=3, data=congo)
summary(m.liar)

m.liar.2 <- ictreg(Y ~ 1, treat="treatment", J=3, data=congo, floor=TRUE)
summary(m.liar.2)

share <- invlogit(m.liar.2$par.treat)-invlogit(m.liar$par.treat)

sim.coefs.congo <- rep(NA, 1000)
#sim.ses.congo <- rep(NA, 1000)

set.seed(140823)
for(i in 1:1000){
  congo$liar <-  rbinom(1000, 1, share)
  congo$wsv.sim <- ifelse(congo$liar==1 & congo$treatment==1 & congo$Y!=4, congo$Y+1, congo$Y)
  
  try(m.out.1.a <- ictreg.joint(wsv.sim ~ female + age.z + edu.z + income.z  + hh_size.z +  murder_yes 
                                + terr_2 + terr_3 + terr_4 + terr_5 + terr_6 + activity_prev,  
                                treat="treatment", 
                                outcome="soc_part.2",
                                outcome.reg="logistic",
                                constrained=TRUE,
                                J=3, data=congo), silent=TRUE)
  
  sim.coefs.congo[i]  <- m.out.1.a$par.outcome[14]
  #sim.ses.congo[i]  <- m.out.1.a$se.outcome[14]
  print(i)
}

#saveRDS(sim.coefs.congo, "sim.coefs.congo.rds")
#saveRDS(sim.ses.congo, "sim.ses.congo.rds")

mean(sim.coefs.congo)
round(quantile(sim.coefs.congo, c(.025, .975)), 2)

par(mar=c(3,3,1,1), mfrow=c(3,1), mgp=c(1.5, .5, 0), cex.axis=.8, yaxs="i", tck=-.005)
hist(sim.coefs.congo, ylim=c(0,350), xlim=c(0,2), main="DRC", ylab="Frequency", xlab="Simulated Effects (Logits) of Wartime Sexual Violence on Civic Engagement", border=F)
#grid(lty=1, lwd=.5, col="lightgrey")
box(col="darkgrey")
abline(v=c(1.11), lty=2, col="red")
abline(v=mean(sim.coefs.congo), lty=2)
text(1.11, 180, "Original Estimate: 1.11", col="red", cex=.8, pos=4)
text(1.065495, 180, "Average of Simulated Estimates: 1.07",  cex=.8, pos=2)

#########################################################################################
### Liberia

m.liar <- ictreg(Y ~ 1, treat="treatment", J=3, data=liberia)
summary(m.liar)

m.liar.2 <- ictreg(Y ~ 1, treat="treatment", J=3, data=liberia, floor=TRUE)
summary(m.liar.2)

share <- invlogit(m.liar.2$par.treat)-invlogit(m.liar$par.treat)

sim.coefs.liberia <- rep(NA, 1000)
#sim.ses.liberia <- rep(NA, 1000)

set.seed(140823)
for(i in 1:1000){
  liberia$liar <-  rbinom(7854, 1, .036) # use same rate as in Congo
  liberia$wsv.sim <- ifelse(liberia$liar==1 & liberia$treatment==1 & liberia$Y!=4, liberia$Y+1, liberia$Y)
  
  try(m.out.1.a <- ictreg.joint(wsv.sim ~ female + age.z + edu.z + income.z + hh_size.z + cw_kill +  as.factor(county), 
                                treat="treatment", 
                                outcome="outcome_ca",
                                outcome.reg="logistic",
                                constrained=TRUE,
                                J=3, data=liberia), silent=TRUE)
  
  sim.coefs.liberia[i]  <- m.out.1.a$par.outcome[10]
  #sim.ses.liberia[i]  <- m.out.1.a$se.outcome[10]
  print(i)
}

#saveRDS(sim.coefs.liberia, "sim.coefs.liberia.rds")
#saveRDS(sim.ses.liberia, "sim.ses.liberia.rds")

mean(sim.coefs.liberia)
round(quantile(sim.coefs.liberia, c(.025, .975)), 2)

par(mar=c(3,3,1,1), mgp=c(1.5, .5, 0), cex.lab=.8, cex.axis=.8, yaxs="i", tck=-.005, las=1)
hist(sim.coefs.liberia, ylim=c(0,350), xlim=c(-.5,2), main="Liberia", ylab="Frequency", xlab="Simulated Effects (Logits) of Wartime Sexual Violence on Civic Engagement", border=F)
#grid(lty=1, lwd=.5, col="lightgrey")
box(col="darkgrey")
abline(v=c(1.75), lty=2, col="red")
abline(v=mean(sim.coefs.liberia), lty=2)
text(1.75, 180, "Original Estimate: 1.75", col="red", cex=.8, pos=2)
text( 0.8018916, 180, "Average of Simulated Estimates: 0.80",  cex=.8, pos=2)

#########################################################################################
### Sri Lanka

m.liar <- ictreg(Y ~ 1, treat="treatment", J=3, data=sri)
summary(m.liar)

m.liar.2 <- ictreg(Y ~ 1, treat="treatment", J=3, data=sri, floor=TRUE)
summary(m.liar.2)

share <- invlogit(m.liar.2$par.treat)-invlogit(m.liar$par.treat)

sim.coefs.sri <- rep(NA, 1000)
#sim.ses.sri <- rep(NA, 1000)

set.seed(140823)
for(i in 1:1000){
  sri$liar <-  rbinom(1800, 1, share) 
  sri$wsv.sim <- ifelse(sri$liar==1 & sri$treatment==1 & sri$Y!=4, sri$Y+1, sri$Y)
  
  try(m.out.1.a <- ictreg.joint(wsv.sim ~ female + age.z + edu.z + income.z + hh_size.z +  killed + trauma + as.factor(Province) + prior,
                                treat="treatment", 
                                outcome="soc_part.2",
                                outcome.reg="logistic",
                                constrained=TRUE,
                                J=3, data=sri), silent=TRUE)
  
  sim.coefs.sri[i]  <- m.out.1.a$par.outcome[18]
  #sim.ses.sri[i]  <- m.out.1.a$se.outcome[18]
  print(i)
}

#saveRDS(sim.coefs.sri, "sim.coefs.sri.rds")
#saveRDS(sim.ses.sri, "sim.ses.sri.rds")

mean(sim.coefs.sri)
round(quantile(sim.coefs.sri, c(.025, .975)), 2)

par(mar=c(3,3,1,1), mgp=c(1.5, .5, 0), cex.lab=.8, cex.axis=.8, yaxs="i", tck=-.005, las=1)
hist(sim.coefs.sri, ylim=c(0,350), xlim=c(-.5,2), main="Sri Lanka", ylab="Frequency", xlab="Simulated Effects (Logits) of Wartime Sexual Violence on Civic Engagement", border=F)
#grid(lty=1, lwd=.5, col="lightgrey")
box(col="darkgrey")
abline(v=c(1.5), lty=2, col="red")
abline(v=mean(sim.coefs.sri), lty=2)
text(1.5, 180, "Original Estimate: 1.50", col="red", cex=.8, pos=2)
text(mean(sim.coefs.sri), 180, "Average of Simulated Estimates: 0.71",  cex=.8, pos=2)

#########################################################################################




